One common way to characterize the structure of a universe is to measure the statistical properties of its haloes. The simplest statistic we can perform is to count the number of haloes found for a given mass interval, normalized to the volume of the universe that contains those structures. This function is usually called the “halo mass function”. In the first part, you will compute the halo mass function for three N-body simulations (DM-only) run using a Λ-CDM cosmological scenario. You will contrast the outputs with theoretical predictions and models. For this activity you will need to install colossus using the command line: pip install colossus
This was uploaded to the Drive of the course and downloaded the full directory (NumEx1), the directory contains
import h5py
import matplotlib as mpl
import matplotlib.pyplot as plt
import matplotlib.gridspec as gridspec
import numpy as np
import pandas as pd
from scipy.interpolate import interp1d
import matplotlib.ticker as ticker
%matplotlib inline
COLOR = 'lightgray'
mpl.rcParams['text.color'] = COLOR # labels, titles, etc colors
mpl.rcParams['axes.labelcolor'] = COLOR # Numbers in axis colors
mpl.rcParams['xtick.color'] = COLOR # ticks colors
mpl.rcParams['ytick.color'] = COLOR # ticks colors
mpl.rc( 'axes',edgecolor = COLOR) # axis Colors
mpl.rc('font', size=15) # Font Size
mpl.rcParams['figure.figsize'] = [12, 12] # Size of figures
# Function to get the values for the extend in plt.imshow
def EXTENT(x,y):
w1=(x[1]-x[0])/2.0
w2=(y[1]-y[0])/2.0
return (x.min()-w1,x.max()+w1,y.min()-w2,y.max()+w2)
#More info, see: https://www.cosmosim.org/metadata/mdr1/rockstar/
file_path = 'mdr1_rockstar.csv'
df = pd.read_csv(file_path)
h = 0.6777
MDR_Mvir = df[['mvir']]/h #Msol/h
MDR_X = df[['x']] #Mpc/h
MDR_Y = df[['y']] #Mpc/h
MDR_Z = df[['z']] #Mpc/h
Limit = (MDR_Mvir>1e10)
Volume= (1000/h)**3
df_filter = df[(df['z'] > 300.) & (df['z'] < 310.)]
Plot of a slice of the dark matter haloes
#plt.hist2d(df_filter['x'],df_filter['y'],bins=200,cmap='hot')
H,x,y=np.histogram2d(df_filter['x'],df_filter['y'],bins=200)
IM=plt.imshow(np.log10(H.T+1),origin='lower',extent=EXTENT(x,y),interpolation='bilinear',cmap='inferno')
CM=plt.colorbar(IM,fraction=0.046, pad=0.04,location='top')
CM.set_label(label='$Log_{10}(Count+1)$',labelpad=20)
plt.axis('scaled')
plt.xlabel('X [cMpc/h]')
plt.ylabel('Y [cMpc/h]')
print()
redshift = [0.0, 0.1, 0.5, 1.0, 2.0]
h = 0.6774
Volume2= (205/h)**3 # Mpc ^3
with h5py.File('TNG300_Dark_z_'+str(redshift[0])+'.hdf5','r') as f1:
TNG300_pos_x = f1['Pos'][:,0]
TNG300_pos_y = f1['Pos'][:,1]
TNG300_pos_z = f1['Pos'][:,2]
TNG300_mvir = f1['M_Crit200'][:]/h
Limit2 = (TNG300_mvir>1e10)
Variables=[]
for i in redshift:
with h5py.File('TNG300_Dark_z_'+str(i)+'.hdf5','r') as f1:
globals()['TNG300_mvir_z_'+str(i)] = f1['M_Crit200'][:]/h
LIM = (globals()['TNG300_mvir_z_'+str(i)]>1e10)
globals()['TNG300_mvir_z_'+str(i)]=globals()['TNG300_mvir_z_'+str(i)][LIM]
Variables.append('TNG300_mvir_z_'+str(i))
Plot of a slice of the dark matter haloes
fil = np.where((TNG300_mvir>1e11)&(TNG300_pos_z>100000)&(TNG300_pos_z<125000))
# plt.hist2d(TNG300_pos_x[fil]/1000.,TNG300_pos_y[fil]/1000.,bins=300,cmap='hot')
H2,x2,y2=np.histogram2d(TNG300_pos_x[fil]/1000.,TNG300_pos_y[fil]/1000.,bins=300)
IM2=plt.imshow((H2.T)**0.75,origin='lower',extent=EXTENT(x2,y2),interpolation='gaussian',cmap='inferno')
CM2=plt.colorbar(IM2,fraction=0.046, pad=0.04,location='top')
CM2.set_label(label='$(Count)^{0.75}$',labelpad=20)
plt.axis('scaled')
plt.xlabel('X [cMpc/h]')
plt.ylabel('Y [cMpc/h]')
print()
file_path = 'EAGLE_100Mpc_DMO.txt'
df_EAGLE = pd.read_csv(file_path, delimiter=',', skiprows=20)
EAGLE_Mvir = df_EAGLE['Group_M_Crit200']#Msol
EAGLE_X = df_EAGLE['GroupCentreOfPotential_x']#Mpc
EAGLE_Y = df_EAGLE['GroupCentreOfPotential_y']#Mpc
EAGLE_Z = df_EAGLE['GroupCentreOfPotential_z']#Mpc
Limit3 = (EAGLE_Mvir>1e10)
Volume3= 100 **3
df_EAGLE_filter = df_EAGLE[(df_EAGLE['GroupCentreOfPotential_z'] > 40.) & (df_EAGLE['GroupCentreOfPotential_z'] < 50.)]
#plt.hist2d(df_EAGLE_filter['GroupCentreOfPotential_x'],df_EAGLE_filter['GroupCentreOfPotential_y'],bins=300,cmap='hot')
H3,x3,y3=np.histogram2d(df_EAGLE_filter['GroupCentreOfPotential_x'],df_EAGLE_filter['GroupCentreOfPotential_y'],bins=300)
IM3=plt.imshow(H3.T,origin='lower',extent=EXTENT(x3,y3),interpolation='bilinear',cmap='inferno')
CM3=plt.colorbar(IM3,fraction=0.046, pad=0.04)
CM3.set_label(label='$Count$',labelpad=20)
plt.axis('scaled')
plt.xlabel('X [cMpc]')
plt.ylabel('Y [cMpc]')
print()
#==================================
# Count particules in 3 dimensions
if True: # High Resolution
x, y, z = EAGLE_X , EAGLE_Y , EAGLE_Z
N= len(x)
#print(N)
H, edges=np.histogramdd((x,y,z),bins=50) #bins=(50,45,55)
x= (edges[0][1:] + edges[0][:-1])/2
y= (edges[1][1:] + edges[1][:-1])/2
z= (edges[2][1:] + edges[2][:-1])/2
x,y,z = np.meshgrid(x,y,z)
x,y,z=x.ravel(),y.ravel(),z.ravel()
#print(len(edges[0]),len(edges[1]),len(edges[2]))
#print(min(x),max(x),min(y),max(y),min(z),max(z))
#print(len(x),len(y),len(z))
if 1: # Scatter density version (plotly)
import plotly.graph_objects as go
from plotly.subplots import make_subplots
#------------------
# Figure
fig = make_subplots(rows=1, cols=1,
specs=[[{'type': 'surface'}]],
print_grid=False,
#subplot_titles=( f'EAGLE'),
)
#------------------
# Scatter plot High resolution
trace1 = go.Volume(name = 'EAGLE',
x = y,
y = x,
z = z,
value = H.ravel(),
#isomin = 0.1,
#isomax = 0.8,
opacity = 0.1, # needs to be small to see through all surfaces
surface_count= 12, # needs to be a large number for good volume rendering
colorscale = 'matter',
colorbar = dict(thickness = 20,
outlinewidth = 0.7,
title = "Count particles per bin",
titleside = "right",
tickmode = "array",
tickvals = [], #[0.001,0.999]
#labelalias = {1: str(H.max()), 0:str(H.min()) },
ticks = "outside",
x = 1.2,) # position of colorbar
)
fig.add_trace(trace1, # trace= : figure
row=1, # row= : Position y
col=1) # col= : Position x
#------------------
# Background color (page and axes)
# https://www.w3schools.com/colors/colors_rgb.asp
fig['layout']['paper_bgcolor']="rgb(17,17,17)"
set_bgcolor= dict(showbackground = True,
backgroundcolor = "rgb(212, 219, 249)", #"rgb(101, 101, 101)"
gridcolor = "rgb(10, 10, 10)",
zeroline = False)
fig.update_scenes(xaxis=set_bgcolor,
yaxis=set_bgcolor,
zaxis=set_bgcolor)
#------------------
# Font color and size
fig.update_layout(font_color = 'white', font_size = 16)
fig.update_layout(height=900, width=900, title_text="EAGLE")
fig.show()
print()
We have defined in the previous step the variables that contain the mass of the halos for each simulation and we also make the conversion to obtain this value without solar masses when applicable.
# Data Counted # This step takes a while
#print(sum(Limit.to_numpy())[0],(MDR_Mvir.to_numpy())[:,0].shape)
h1,x1=np.histogram(np.log10(MDR_Mvir[Limit]),bins=100) # We count the data in 100 bins log spaced
w1= (x1[1]-x1[0]) # Width of the bins (mass enclosed)
X1= (x1[1:]+x1[:-1])/2. # Center of the bins
#print(sum(Limit2),TNG300_mvir.shape)
h2,x2=np.histogram(np.log10(TNG300_mvir[Limit2]),bins=100)
w2= (x2[1]-x2[0])
X2= (x2[1:]+x2[:-1])/2.
#print(sum(Limit3),EAGLE_Mvir.shape)
h3,x3=np.histogram(np.log10(EAGLE_Mvir[Limit3]),bins=100)
w3= (x3[1]-x3[0])
X3= (x3[1:]+x3[:-1])/2.
plt.style.use('dark_background')
plt.style.use('seaborn-bright')
plt.stairs(h1,x1,fill=False,label='Small MultiDark Planck simulation')
plt.stairs(h2,x2,fill=False,label='IllustrisTNG300-Dark')
plt.stairs(h3,x3,fill=False,label='EAGLE DM-ONLY')
plt.yscale('log')
plt.ylabel('$N_{halo}$')
plt.xlabel('$\log_{10}(M [M_\odot])$')
plt.legend()
plt.grid(ls=":",color='lightgray',lw=0.5,which='both')
(The variables of volume were defined when we read the data)
plt.plot(X1,h1/w1/Volume,marker='o',label='Small MultiDark Planck simulation')
plt.plot(X2,h2/w2/Volume2,marker='d',label='IllustrisTNG300-Dark')
plt.plot(X3,h3/w3/Volume3,marker='s',label='EAGLE DM-ONLY')
plt.yscale('log')
#plt.xscale('log')
plt.ylabel(r'$N_{halo}\cdot\ dv^{-1}\cdot\ d(Log_{10}(M_{200}))^{-1}$',fontsize=20)
plt.xlabel('$\log_{10}(M [M_\odot])$')
plt.legend(loc='upper left', bbox_to_anchor=(1.05, 1.))
plt.grid(ls=":",color='lightgray',lw=0.5,which='both')
from colossus.cosmology import cosmology
from colossus.lss import mass_function
we will create the distributions and make an interpolation of them in order to make much easier the comparison
cosmology.setCosmology('WMAP9')
new_x = np.linspace(9.99,16.01,100)
mfunc_P = mass_function.massFunction(10**new_x, 0.0, mdef = 'fof', model = 'press74',q_out='dndlnM')
mfunc_S = mass_function.massFunction(10**new_x, 0.0, mdef = 'fof', model = 'sheth99',q_out='dndlnM')
# EXTRA = mass_function.massFunction(10**new_x, 0.0, mdef = 'vir', model = 'tinker08',q_out='dndlnM')
MH1= interp1d(new_x,mfunc_P)
MH2= interp1d(new_x,mfunc_S)
Y1=h1/w1/Volume
Y1[Y1<=0]=np.nan # easier to handle
Y2=h2/w2/Volume2
Y2[Y2<=0]=np.nan # easier to handle
Y3=h3/w3/Volume3
Y3[Y3<=0]=np.nan # easier to handle
plt.plot(X1,Y1,lw=3,label='Small MultiDark Planck simulation',color='b')#,marker='o')
plt.plot(X2,Y2,lw=3,label='IllustrisTNG300-Dark',color='springgreen')#,marker='d')
plt.plot(X3,Y3,lw=3,label='EAGLE DM-ONLY',color='r')#,marker='s')
plt.plot(new_x,mfunc_P,ls='--',lw=3,label='Press and Schechter',color='fuchsia')#,marker='.')
plt.plot(new_x,mfunc_S,ls='--',lw=3,label='Sheth and Tormen',color='aqua')#,marker='*')
# plt.plot(new_x,EXTRA,label='EXTRA',marker='p')
plt.ylabel(r'$N_{halo}\cdot\ dv^{-1}\cdot\ d(Log_{10}(M_{200}))^{-1}$',fontsize=20)
plt.xlabel('$\log_{10}(M [M_\odot])$')
plt.legend(loc='upper left', bbox_to_anchor=(1.02, 1.))
plt.grid(ls=":",color='lightgray',lw=0.5,which='both')
plt.ylim(2e-8,2)
plt.yscale('log')
fig= plt.figure(figsize=(18,9))
gs = gridspec.GridSpec( 3,3, # Number of axis y,x
height_ratios = [0.5,0.05,0.5], # relatives ratios of heigh
width_ratios = [0.33,0.33,0.33], # relatives ratio of with
left = 0.1, # Space to the edge of left from the nearest axis
right = 0.95, # Space to the edge of right from the nearest axis
bottom= 0.1, # Space to the edge of bottom from the nearest axis
top = 0.97, # Space to the edge of top from the nearest axis
wspace= 0.15, # Space horizontal between each of the axis
hspace= 0.0) # Space vertical between each of the axis
ax1 = fig.add_subplot(gs[0,0])
ax2 = fig.add_subplot(gs[2,0],sharex=ax1,sharey=ax1)
ax3 = fig.add_subplot(gs[0,1],sharex=ax1,sharey=ax1)
ax4 = fig.add_subplot(gs[2,1],sharex=ax1,sharey=ax1)
ax5 = fig.add_subplot(gs[0,2],sharex=ax1,sharey=ax1)
ax6 = fig.add_subplot(gs[2,2],sharex=ax1,sharey=ax1)
# (here we can define if we want to share axes and with which one)
#--
# Ploting
Center=1.0
Op= lambda a,b: a/b
ax1.set_yscale('log')
ax1.yaxis.set_major_formatter(ticker.FuncFormatter(lambda y,pos: ('{{:.{:1d}f}}'.format(int(np.maximum(-np.log10(y),0)))).format(y)))
fig.supxlabel('$\log_{10}(M [M_\odot])$',fontsize=18)
fig.supylabel('Ratio',fontsize=18)
#fig.supylabel('Residual',fontsize=18)
# Press and Schechter
ax1.plot(X1,Op(Y1,MH1(X1)),lw=3,label='Small MultiDark Planck simulation',color='b')
ax1.axhline(Center,ls='--',lw=3,label='Press and Schechter',color='fuchsia',alpha=0.5)
ax1.plot([],[],ls='--',lw=3,label='Sheth and Tormen',color='aqua',alpha=0.5)
ax1.legend(loc='upper center', bbox_to_anchor=(.5, 1.25))
ax3.plot(X2,Op(Y2,MH1(X2)),lw=3,label='IllustrisTNG300-Dark',color='springgreen')
ax3.axhline(Center,ls='--',lw=3,label='Press and Schechter',color='fuchsia',alpha=0.5)
ax3.plot([],[],ls='--',lw=3,label='Sheth and Tormen',color='aqua',alpha=0.5)
ax3.legend(loc='upper center', bbox_to_anchor=(.5, 1.25))
ax5.plot(X3,Op(Y3,MH1(X3)),lw=3,label='EAGLE DM-ONLY',color='r')
ax5.axhline(Center,ls='--',lw=3,label='Press and Schechter',color='fuchsia',alpha=0.5)
ax5.plot([],[],ls='--',lw=3,label='Sheth and Tormen',color='aqua',alpha=0.5)
ax5.legend(loc='upper center', bbox_to_anchor=(.5, 1.25))
# Sheth and Tormen
ax2.plot(X1,Op(Y1,MH2(X1)),lw=3,label='Small MultiDark Planck simulation',color='b')
ax2.axhline(Center,ls='--',lw=3,label='Sheth and Tormen',color='aqua',alpha=0.5)
# ax2.legend(loc='upper right')
ax4.plot(X2,Op(Y2,MH2(X2)),lw=3,label='IllustrisTNG300-Dark',color='springgreen')
ax4.axhline(Center,ls='--',lw=3,label='Sheth and Tormen',color='aqua',alpha=0.5)
# ax4.legend(loc='upper right')
ax6.plot(X3,Op(Y3,MH2(X3)),lw=3,label='EAGLE DM-ONLY',color='r')
ax6.axhline(Center,ls='--',lw=3,label='Sheth and Tormen',color='aqua',alpha=0.5)
# ax6.legend(loc='upper right')
# Erasing extra y labels
plt.setp(ax1.get_xticklabels(), visible=False)
plt.setp(ax3.get_xticklabels(), visible=False)
plt.setp(ax5.get_xticklabels(), visible=False)
print()
From this last plot we can clearly see that the ellipsoidal collapse fits better over all the simulation, now if we look more in detail:
# We read the variables (they were defined when we read the TNG-300-dark-1)
Sec_Var={}
for i in Variables:
D=globals()[i]
#print(D.shape)
HH,XX=np.histogram(np.log10(D),bins=100)
mask= (HH>0)
WW = (XX[1]-XX[0])
XX = ((XX[1:]+XX[:-1])/2.)[mask]
Sec_Var[i]=np.array([XX,HH[mask]/WW/Volume2])
fig= plt.figure(figsize=(18,9))
gs = gridspec.GridSpec( 1,3, # Number of axis y,x
height_ratios = [1], # relatives ratios of heigh
width_ratios = [0.5,0.1,0.5], # relatives ratio of with
left = 0.1, # Space to the edge of left from the nearest axis
right = 0.95, # Space to the edge of right from the nearest axis
bottom= 0.1, # Space to the edge of bottom from the nearest axis
top = 0.97, # Space to the edge of top from the nearest axis
wspace= 0.15, # Space horizontal between each of the axis
hspace= 0.0) # Space vertical between each of the axis
ax1 = fig.add_subplot(gs[0,0])
ax2 = fig.add_subplot(gs[0,2],sharex=ax1,sharey=ax1)
print('File\t\t\t[Slope (N) Coeff (M0)]')
print('-------------------------------------------------')
for i in Sec_Var.keys():
D = Sec_Var[i]
ax1.plot(D[0],D[1])
pf=np.polyfit(D[0],np.log10(D[1]),1)
print(f'{i}\t{pf}')
ax2.plot(D[0],10**((D[0]*pf[0])+pf[1]),label=i)
ax1.set_yscale('log')
ax2.legend()
ax1.grid(zorder=0,ls=':',color='gray')
ax2.grid(zorder=0,ls=':',color='gray')
fig.supylabel(r'$N_{halo}\cdot\ dv^{-1}\cdot\ d(Log_{10}(M_{200}))^{-1}$',fontsize=20)
fig.supxlabel('$\log_{10}(M [M_\odot])$',fontsize=18)
print()
File [Slope (N) Coeff (M0)] ------------------------------------------------- TNG300_mvir_z_0.0 [-1.01584723 9.7721955 ] TNG300_mvir_z_0.1 [-1.01394106 9.75779471] TNG300_mvir_z_0.5 [-1.06472794 10.3614411 ] TNG300_mvir_z_1.0 [-1.12013746 10.98513452] TNG300_mvir_z_2.0 [-1.24970781 12.3532928 ]
The slope presents slighty changues over redshift , in specific this decrease with redshift and this makes sense if we assume the merger trees evolution of the universe (we have less massive halos at early stages of the universe and more of small ones)
masses = [1e10, 1e11, 1e12, 1e13]
Data = []
for d in Variables: # Getting the info of the specific redshift
D = globals()[d]
aux = []
for i in range(len(masses)-1):
l,u = masses[i],masses[i+1] # Range of masses desired
mask = (l<D)&(D<u)
aux.append(len(D[mask])/Volume2/(np.log10(u)-np.log10(l))) # Number of halos/ Volume / bin mass width
Data.append(aux)
Data=np.array(Data)
fig= plt.figure(figsize=(18,6))
gs = gridspec.GridSpec( 1,5, # Number of axis y,x
height_ratios = [1], # relatives ratios of heigh
width_ratios = [1,0.1,1,0.1,1], # relatives ratio of with
left = 0.1, # Space to the edge of left from the nearest axis
right = 0.95, # Space to the edge of right from the nearest axis
bottom= 0.15, # Space to the edge of bottom from the nearest axis
top = 0.97, # Space to the edge of top from the nearest axis
wspace= 0.1, # Space horizontal between each of the axis
hspace= 0.0) # Space vertical between each of the axis
ax1 = fig.add_subplot(gs[0,0])
ax2 = fig.add_subplot(gs[0,2])
ax3 = fig.add_subplot(gs[0,4])
AX=[ax1,ax2,ax3]
# Function to put scientific notation
def nsf(x):
return f'10^{{{np.log10(x):.0f}}}'
# Ploting the masses
for i in range(0,len(masses)-1):
l,u = masses[i],masses[i+1] # Range of masses desired
DDD = Data[:,i]
AX[i].plot(redshift,DDD,marker='o',label=f'${nsf(l)}-{nsf(u)}$')
AX[i].legend(loc='upper left', bbox_to_anchor=(0.52, 1.12))
AX[i].ticklabel_format(axis='y', scilimits=[-1, 1])
AX[i].grid(zorder=0,ls=':',color='gray')
#AX[i].set_yscale('log')
fig.supylabel(r'$N_{halo}\cdot\ dv^{-1}\cdot\ d(Log_{10}(M_{200}))^{-1}$',fontsize=20)
fig.supxlabel('Redshift',fontsize=20)
print()
This behaviour can be easily understood in the context of the hierarchical universe, this means that the amount of small halos (first panel) will naturally decreases as function of redshift due they merge along the time , in the other extreme the most massive ones (last panel) will be more abundant in later stages of the universe since more halos have had enough time to being formed, finally the middle panel reach a peak at certain redshift where they have just enough time to being formed but not much to be merged.
print('Datapoint in the last graph')
Data
Datapoint in the last graph
array([[0.10847032, 0.01450199, 0.001861 ],
[0.11198504, 0.0149035 , 0.00189412],
[0.12394938, 0.0161267 , 0.00197 ],
[0.13284127, 0.0167439 , 0.00188525],
[0.13456275, 0.01507434, 0.00128826]])